Code for calculating climate distance (Figure S3; Table S4)

To Do

Load Libraries

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(boot) #for bootstrapping 
library(broom) #for tidy function
library(geosphere) #for calculating geographic distance

Load Climate Data

all_clim <- read_csv("../Processed.Data/Climate/All_Clim.csv")
## Rows: 2951 Columns: 23
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (4): parent.pop, elevation.group, timeframe, Season
## dbl (19): elev_m, Lat, Long, year, cwd, pck, ppt, tmn, tmx, ann_tmean, mean_...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

WL2 Garden Climate

#Nov 2022-October 2023
wl2_garden_2023_wtryr <- all_clim %>% 
  filter(parent.pop=="WL2_Garden", Season=="Water Year", year==2023) %>%
  rename_with(.cols = c(8:22), ~paste0(.x, '_WL2')) %>% 
  select(cwd_WL2:ppt_coldest_month_WL2)
#Nov 2023-October 2024
wl2_garden_2024_wtryr <- all_clim %>% 
  filter(parent.pop=="WL2_Garden", Season=="Water Year", year==2024) %>%
  rename_with(.cols = c(8:22), ~paste0(.x, '_WL2')) %>% 
  select(cwd_WL2:ppt_coldest_month_WL2)
#July-December 2023 
wl2_garden_2023_grwssn <- all_clim %>% 
  filter(parent.pop=="WL2_Garden", Season=="Growth Season", year==2023) %>%
  rename_with(.cols = c(8:22), ~paste0(.x, '_WL2')) %>% 
  select(cwd_WL2:ppt_coldest_month_WL2)
#June-November 2024
wl2_garden_2024_grwssn <- all_clim %>% 
  filter(parent.pop=="WL2_Garden", Season=="Growth Season", year==2024) %>%
  rename_with(.cols = c(8:22), ~paste0(.x, '_WL2')) %>% 
  select(cwd_WL2:ppt_coldest_month_WL2)

Gower’s Climate Distance

(1/P) * SUM ((absolute value(Ai - Bi)) / range(i)) for each variable

  • P = number of environmental variables = 13 (without CWD)

  • Ai = 30 year avg of that variable for the home site

  • Bi = July 2023-Dec 2023 avg of that variable for the WL2 garden

  • Range(i) = maximum - minimum of that variable in the whole data set (across sites)

Prep the dataframes

#pop info 
pops <- all_clim %>% 
  filter(!is.na(Lat), parent.pop!="WL2_Garden") %>% 
  select(parent.pop, elevation.group, elev_m, Lat, Long) %>% 
  unique()

#home climates
recent_wtryr_clim <- all_clim %>% 
  filter(timeframe=="recent", Season=="Water Year")
recent_wtryr_clim_nest <- recent_wtryr_clim %>% nest(.by=year) #nest to prepare for bootstrapping 
recent_wtryr_clim_nest
## # A tibble: 30 × 2
##     year data              
##    <dbl> <list>            
##  1  1994 <tibble [23 × 22]>
##  2  1995 <tibble [23 × 22]>
##  3  1996 <tibble [23 × 22]>
##  4  1997 <tibble [23 × 22]>
##  5  1998 <tibble [23 × 22]>
##  6  1999 <tibble [23 × 22]>
##  7  2000 <tibble [23 × 22]>
##  8  2001 <tibble [23 × 22]>
##  9  2002 <tibble [23 × 22]>
## 10  2003 <tibble [23 × 22]>
## # ℹ 20 more rows
historic_wtryr_clim <- all_clim %>% 
  filter(timeframe=="historic", Season=="Water Year")
historic_wtryr_clim_nest <- historic_wtryr_clim %>% nest(.by=year) #nest to prepare for bootstrapping 
historic_wtryr_clim_nest
## # A tibble: 30 × 2
##     year data              
##    <dbl> <list>            
##  1  1964 <tibble [23 × 22]>
##  2  1965 <tibble [23 × 22]>
##  3  1966 <tibble [23 × 22]>
##  4  1967 <tibble [23 × 22]>
##  5  1968 <tibble [23 × 22]>
##  6  1969 <tibble [23 × 22]>
##  7  1970 <tibble [23 × 22]>
##  8  1971 <tibble [23 × 22]>
##  9  1972 <tibble [23 × 22]>
## 10  1973 <tibble [23 × 22]>
## # ℹ 20 more rows
recent_grwssn_clim <- all_clim %>% 
  filter(timeframe=="recent", Season=="Growth Season", parent.pop!="WL2_Garden") 
recent_grwssn_clim_nest <- recent_grwssn_clim %>% nest(.by=year) #nest to prepare for bootstrapping 
recent_grwssn_clim_nest
## # A tibble: 30 × 2
##     year data              
##    <dbl> <list>            
##  1  1994 <tibble [23 × 22]>
##  2  1995 <tibble [23 × 22]>
##  3  1996 <tibble [23 × 22]>
##  4  1997 <tibble [23 × 22]>
##  5  1998 <tibble [23 × 22]>
##  6  1999 <tibble [23 × 22]>
##  7  2000 <tibble [23 × 22]>
##  8  2001 <tibble [23 × 22]>
##  9  2002 <tibble [23 × 22]>
## 10  2003 <tibble [23 × 22]>
## # ℹ 20 more rows
historic_grwssn_clim <- all_clim %>% 
  filter(timeframe=="historic", Season=="Growth Season", parent.pop!="WL2_Garden")
historic_grwssn_clim_nest <- historic_grwssn_clim %>% nest(.by=year) #nest to prepare for bootstrapping 
historic_grwssn_clim_nest
## # A tibble: 30 × 2
##     year data              
##    <dbl> <list>            
##  1  1964 <tibble [23 × 22]>
##  2  1965 <tibble [23 × 22]>
##  3  1966 <tibble [23 × 22]>
##  4  1967 <tibble [23 × 22]>
##  5  1968 <tibble [23 × 22]>
##  6  1969 <tibble [23 × 22]>
##  7  1970 <tibble [23 × 22]>
##  8  1971 <tibble [23 × 22]>
##  9  1972 <tibble [23 × 22]>
## 10  1973 <tibble [23 × 22]>
## # ℹ 20 more rows
#WL2 range prep
wl2_garden_2023_wtryr_range_prep <- wl2_garden_2023_wtryr %>% 
  rename_with(~str_remove(., "_WL2"), everything())
wl2_garden_2024_wtryr_range_prep <- wl2_garden_2024_wtryr %>% 
  rename_with(~str_remove(., "_WL2"), everything())
wl2_garden_2023_grwssn_range_prep <- wl2_garden_2023_grwssn %>% 
  rename_with(~str_remove(., "_WL2"), everything())
wl2_garden_2024_grwssn_range_prep <- wl2_garden_2024_grwssn %>% 
  rename_with(~str_remove(., "_WL2"), everything())

Create the gower_calc function

#gower's calc function 
gowers_calc <- function(data, indices, P, garden_clim, garden_range) { #function with all of the code necessary for calculating gowers distance 
  #data = _clim_nest (recent or historical) - needs to be nested by year; P = # climate variables;
  #garden_clim = wl2 garden clim with _Wl2 at end of column names
  #garden_range = wl2 garden data with WL2 removed from the column names 
  
  data <-data[indices,] # subset per bootstrap indices
  
  data <- data %>% unnest(data) #unnest so the function can access the climate data
  
  data_means <- data %>% 
    group_by(parent.pop, elevation.group, elev_m, Lat, Long, timeframe) %>% 
    summarise_at(c("cwd",  "pck", "ppt", "tmn", "tmx", "ann_tmean", "mean_diurnal_range", 
                   "temp_seasonality", "temp_ann_range",
                 "tmean_wettest_month", "tmean_driest_month", "ann_ppt",
                 "ppt_seasonality","ppt_warmest_month", "ppt_coldest_month"),
               c(mean), na.rm = TRUE) #get 30 year averages for each climate variable 
  
  range_merge <- bind_rows(data_means, garden_range)
  
  WL2_home_climate_ranges <- range_merge %>% #calculate ranges
    ungroup() %>% 
  summarise(cwd_range=max(cwd)-min(cwd),
            pck_range=max(pck)-min(pck),
            ppt_range=max(ppt)-min(ppt), 
            tmn_range=max(tmn)-min(tmn), 
            tmx_range=max(tmx)-min(tmx), 
            ann_tmean_range=max(ann_tmean)-min(ann_tmean),
            mean_diurnal_range_range=max(mean_diurnal_range)-min(mean_diurnal_range),
            temp_seasonality_range=max(temp_seasonality)-min(temp_seasonality),
            temp_ann_range_range=max(temp_ann_range)-min(temp_ann_range),
            tmean_wettest_month_range=max(tmean_wettest_month)-min(tmean_wettest_month),
            tmean_driest_month_range=max(tmean_driest_month)-min(tmean_driest_month),
            ann_ppt_range=max(ann_ppt)-min(ann_ppt), 
            ppt_seasonality_range=max(ppt_seasonality)-min(ppt_seasonality),
            ppt_warmest_month_range=max(ppt_warmest_month)-min(ppt_warmest_month), 
            ppt_coldest_month_range=max(ppt_coldest_month)-min(ppt_coldest_month))
  
  WL2_home_climate <- bind_cols(garden_clim, data_means) #add WL2 climate data to home climate data 
  
  WL2_home_climate_with_ranges <- bind_cols(WL2_home_climate, WL2_home_climate_ranges) #add in ranges 
  
  gowers_calc_each_var <- WL2_home_climate_with_ranges %>% #variable by variable calc
  mutate(cwd_gowers=abs(cwd_WL2-cwd) / cwd_range,
         pck_gowers=abs(pck_WL2-pck) / pck_range,
         ppt_gowers=abs(ppt_WL2 - ppt) / ppt_range,
         tmn_gowers=abs(tmn_WL2 - tmn) / tmn_range,
         tmx_gowers=abs(tmx_WL2 - tmx) / tmx_range,
         ann_tmean_gowers=abs(ann_tmean_WL2 - ann_tmean) / ann_tmean_range,
         mean_diurnal_range_gowers=abs(mean_diurnal_range_WL2 - mean_diurnal_range) / mean_diurnal_range_range,
         temp_seasonality_gowers=abs(temp_seasonality_WL2 - temp_seasonality) / temp_seasonality_range,
         temp_ann_range_gowers=abs(temp_ann_range_WL2 - temp_ann_range) / temp_ann_range_range,
         tmean_wettest_month_gowers=abs(tmean_wettest_month_WL2 - tmean_wettest_month) / tmean_wettest_month_range,
         tmean_driest_month_gowers=abs(tmean_driest_month_WL2 - tmean_driest_month) / tmean_driest_month_range,
         ann_ppt_gowers=abs(ann_ppt_WL2 - ann_ppt) / ann_ppt_range,
         ppt_seasonality_gowers=abs(ppt_seasonality_WL2 - ppt_seasonality) / ppt_seasonality_range,
         ppt_warmest_month_gowers=abs(ppt_warmest_month_WL2 - ppt_warmest_month) / ppt_warmest_month_range,
         ppt_coldest_month_gowers=abs(ppt_coldest_month_WL2 - ppt_coldest_month) / ppt_coldest_month_range) %>% 
  dplyr::select(parent.pop, elevation.group, elev_m, ends_with("_gowers"))

 gowers_calc_per_pop <- gowers_calc_each_var %>% #final gowers calc 
  mutate(Gowers_Dist=(1/P)*(cwd_gowers + pck_gowers + ppt_gowers + tmn_gowers + tmx_gowers +
                                ann_tmean_gowers + mean_diurnal_range_gowers +
                                temp_seasonality_gowers +temp_ann_range_gowers +
                                tmean_wettest_month_gowers +
                                tmean_driest_month_gowers +ann_ppt_gowers +
                                ppt_seasonality_gowers + ppt_warmest_month_gowers +
                                ppt_coldest_month_gowers)) %>% 
  dplyr::select(parent.pop, elevation.group, elev_m, Gowers_Dist)
  
 gowers_calc_per_pop %>% pull(Gowers_Dist) #make the result a vector 
}

Bootstrapping

2023 Garden

# water year - recent
gowers.boot_2023_recent_wtryr <- boot(data=recent_wtryr_clim_nest, statistic=gowers_calc, R=1000, P=15, garden_clim=wl2_garden_2023_wtryr, garden_range=wl2_garden_2023_wtryr_range_prep) #will sample each row (year) with replacement 

for(i in 1:23) {
  plot(gowers.boot_2023_recent_wtryr, index=i) #distributions look normal for the most part 
}

boot_recent_wtryr_results_2023 <- tidy(gowers.boot_2023_recent_wtryr,conf.int=TRUE,conf.method="norm") %>%  #all pops
  rename(Gowers_Dist = statistic) %>% 
  mutate(timeframe="Recent", Season="Water Year", Year=2023)
boot_gowers_recent_wtryr_pops_2023 <- bind_cols(pops, boot_recent_wtryr_results_2023) %>% 
  select(parent.pop:Long, timeframe, Season, Year, Gowers_Dist, conf.low, conf.high) %>% 
  arrange(Gowers_Dist)

# water year - historic
gowers.boot_2023_historic_wtryr <- boot(data=historic_wtryr_clim_nest, statistic=gowers_calc, R=1000, P=15, garden_clim=wl2_garden_2023_wtryr, garden_range=wl2_garden_2023_wtryr_range_prep) #will sample each row (year) with replacement 

for(i in 1:23) {
  plot(gowers.boot_2023_historic_wtryr, index=i) #distributions look normal for the most part 
}

boot_historic_wtryr_results_2023 <- tidy(gowers.boot_2023_historic_wtryr,conf.int=TRUE,conf.method="norm") %>%  #all pops
  rename(Gowers_Dist = statistic) %>% 
  mutate(timeframe="Historic", Season="Water Year", Year=2023)
boot_gowers_historic_wtryr_pops_2023 <- bind_cols(pops, boot_historic_wtryr_results_2023) %>% 
  select(parent.pop:Long, timeframe, Season, Year, Gowers_Dist, conf.low, conf.high) %>% 
  arrange(Gowers_Dist)

# growth season - recent
gowers.boot_2023_recent_grwssn <- boot(data=recent_grwssn_clim_nest, statistic=gowers_calc, R=1000, P=15, garden_clim=wl2_garden_2023_grwssn, garden_range=wl2_garden_2023_grwssn_range_prep) #will sample each row (year) with replacement 

for(i in 1:23) {
  plot(gowers.boot_2023_recent_grwssn, index=i) #distributions look normal for the most part 
}

boot_recent_grwssn_results_2023 <- tidy(gowers.boot_2023_recent_grwssn,conf.int=TRUE,conf.method="norm") %>%  #all pops
  rename(Gowers_Dist = statistic) %>% 
  mutate(timeframe="Recent", Season="Growth Season", Year=2023)
boot_gowers_recent_grwssn_pops_2023 <- bind_cols(pops, boot_recent_grwssn_results_2023) %>% 
  select(parent.pop:Long, timeframe, Season, Year, Gowers_Dist, conf.low, conf.high) %>% 
  arrange(Gowers_Dist)

# growth season - historic
gowers.boot_2023_historic_grwssn <- boot(data=historic_grwssn_clim_nest, statistic=gowers_calc, R=1000, P=15, garden_clim=wl2_garden_2023_grwssn, garden_range=wl2_garden_2023_grwssn_range_prep) #will sample each row (year) with replacement 

for(i in 1:23) {
  plot(gowers.boot_2023_historic_grwssn, index=i) #distributions look normal for the most part 
}

boot_historic_grwssn_results_2023 <- tidy(gowers.boot_2023_historic_grwssn,conf.int=TRUE,conf.method="norm") %>%  #all pops
  rename(Gowers_Dist = statistic) %>% 
  mutate(timeframe="Historic", Season="Growth Season", Year=2023)
boot_gowers_historic_grwssn_pops_2023 <- bind_cols(pops, boot_historic_grwssn_results_2023) %>% 
  select(parent.pop:Long, timeframe, Season, Year, Gowers_Dist, conf.low, conf.high) %>% 
  arrange(Gowers_Dist)

2024 Garden

# water year - recent
gowers.boot_2024_recent_wtryr <- boot(data=recent_wtryr_clim_nest, statistic=gowers_calc, R=1000, P=15, garden_clim=wl2_garden_2024_wtryr, garden_range=wl2_garden_2024_wtryr_range_prep) #will sample each row (year) with replacement 

for(i in 1:23) {
  plot(gowers.boot_2024_recent_wtryr, index=i) #distributions look normal for the most part 
}

boot_recent_wtryr_results_2024 <- tidy(gowers.boot_2024_recent_wtryr,conf.int=TRUE,conf.method="norm") %>%  #all pops
  rename(Gowers_Dist = statistic) %>% 
  mutate(timeframe="Recent", Season="Water Year", Year=2024)
boot_gowers_recent_wtryr_pops_2024 <- bind_cols(pops, boot_recent_wtryr_results_2024) %>% 
  select(parent.pop:Long, timeframe, Season, Year, Gowers_Dist, conf.low, conf.high) %>% 
  arrange(Gowers_Dist)

# water year - historic
gowers.boot_2024_historic_wtryr <- boot(data=historic_wtryr_clim_nest, statistic=gowers_calc, R=1000, P=15, garden_clim=wl2_garden_2024_wtryr, garden_range=wl2_garden_2024_wtryr_range_prep) #will sample each row (year) with replacement 

for(i in 1:23) {
  plot(gowers.boot_2024_historic_wtryr, index=i) #distributions look normal for the most part 
}

boot_historic_wtryr_results_2024 <- tidy(gowers.boot_2024_historic_wtryr,conf.int=TRUE,conf.method="norm") %>%  #all pops
  rename(Gowers_Dist = statistic) %>% 
  mutate(timeframe="Historic", Season="Water Year", Year=2024)
boot_gowers_historic_wtryr_pops_2024 <- bind_cols(pops, boot_historic_wtryr_results_2024) %>% 
  select(parent.pop:Long, timeframe, Season, Year, Gowers_Dist, conf.low, conf.high) %>% 
  arrange(Gowers_Dist)

# growth season - recent
gowers.boot_2024_recent_grwssn <- boot(data=recent_grwssn_clim_nest, statistic=gowers_calc, R=1000, P=15, garden_clim=wl2_garden_2024_grwssn, garden_range=wl2_garden_2024_grwssn_range_prep) #will sample each row (year) with replacement 

for(i in 1:23) {
  plot(gowers.boot_2024_recent_grwssn, index=i) #distributions look normal for the most part 
}

boot_recent_grwssn_results_2024 <- tidy(gowers.boot_2024_recent_grwssn,conf.int=TRUE,conf.method="norm") %>%  #all pops
  rename(Gowers_Dist = statistic) %>% 
  mutate(timeframe="Recent", Season="Growth Season", Year=2024)
boot_gowers_recent_grwssn_pops_2024 <- bind_cols(pops, boot_recent_grwssn_results_2024) %>% 
  select(parent.pop:Long, timeframe, Season, Year, Gowers_Dist, conf.low, conf.high) %>% 
  arrange(Gowers_Dist)

# growth season - historic
gowers.boot_2024_historic_grwssn <- boot(data=historic_grwssn_clim_nest, statistic=gowers_calc, R=1000, P=15, garden_clim=wl2_garden_2024_grwssn, garden_range=wl2_garden_2024_grwssn_range_prep) #will sample each row (year) with replacement 

for(i in 1:23) {
  plot(gowers.boot_2024_historic_grwssn, index=i) #distributions look normal for the most part 
}

boot_historic_grwssn_results_2024 <- tidy(gowers.boot_2024_historic_grwssn,conf.int=TRUE,conf.method="norm") %>%  #all pops
  rename(Gowers_Dist = statistic) %>% 
  mutate(timeframe="Historic", Season="Growth Season", Year=2024)
boot_gowers_historic_grwssn_pops_2024 <- bind_cols(pops, boot_historic_grwssn_results_2024) %>% 
  select(parent.pop:Long, timeframe, Season, Year, Gowers_Dist, conf.low, conf.high) %>% 
  arrange(Gowers_Dist)

Average gower’s distance

all_gowers <- bind_rows(boot_gowers_recent_wtryr_pops_2023,
                            boot_gowers_historic_wtryr_pops_2023,
                            boot_gowers_recent_grwssn_pops_2023,
                            boot_gowers_historic_grwssn_pops_2023,
                            boot_gowers_recent_wtryr_pops_2024,
                            boot_gowers_historic_wtryr_pops_2024,
                            boot_gowers_recent_grwssn_pops_2024,
                            boot_gowers_historic_grwssn_pops_2024) 
wl2_gowers_avg <- all_gowers %>% 
  group_by(parent.pop, elevation.group, elev_m, Lat, Long, timeframe, Season) %>% 
  summarise(Avg_Gowers=mean(Gowers_Dist, na.rm=TRUE))
## `summarise()` has grouped output by 'parent.pop', 'elevation.group', 'elev_m',
## 'Lat', 'Long', 'timeframe'. You can override using the `.groups` argument.

Plot the results

# 2023 water year
boot_gowers_recent_wtryr_pops_2023 %>% 
  bind_rows(boot_gowers_historic_wtryr_pops_2023) %>% 
  ggplot(aes(x=fct_reorder(parent.pop, Gowers_Dist), y=Gowers_Dist, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high),
                width=.1, position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="Gowers Envtal Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe)

ggsave("../Figures/wtryear_2023Gowers_fromWL2.png", width = 24, height = 8, units = "in")

# 2023 growth season
boot_gowers_recent_grwssn_pops_2023 %>% 
  bind_rows(boot_gowers_historic_grwssn_pops_2023) %>% 
  ggplot(aes(x=fct_reorder(parent.pop, Gowers_Dist), y=Gowers_Dist, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high),
                width=.1, position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="Gowers Envtal Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe)

ggsave("../Figures/growthseason_2023Gowers_fromWL2.png", width = 24, height = 8, units = "in")

# 2024 water year 
boot_gowers_recent_wtryr_pops_2024 %>% 
  bind_rows(boot_gowers_historic_wtryr_pops_2024) %>% 
  ggplot(aes(x=fct_reorder(parent.pop, Gowers_Dist), y=Gowers_Dist, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high),
                width=.1, position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="Gowers Envtal Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe)

ggsave("../Figures/wtryear_2024Gowers_fromWL2.png", width = 24, height = 8, units = "in")

# 2024 growth season
boot_gowers_recent_grwssn_pops_2024 %>% 
  bind_rows(boot_gowers_historic_grwssn_pops_2024) %>% 
  ggplot(aes(x=fct_reorder(parent.pop, Gowers_Dist), y=Gowers_Dist, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high),
                width=.1, position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="Gowers Envtal Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe)

ggsave("../Figures/growthseason_2024Gowers_fromWL2.png", width = 24, height = 8, units = "in")

# average water year 
wl2_gowers_avg %>% 
  filter(Season=="Water Year") %>% 
  ggplot(aes(x=fct_reorder(parent.pop, Avg_Gowers), y=Avg_Gowers, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="Gowers Envtal Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe)

ggsave("../Figures/wtryear_AvgGowers_fromWL2.png", width = 24, height = 8, units = "in")

# average growth season 
wl2_gowers_avg %>% 
  filter(Season=="Growth Season") %>% 
  ggplot(aes(x=fct_reorder(parent.pop, Avg_Gowers), y=Avg_Gowers, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="Gowers Envtal Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe)

ggsave("../Figures/growthseason_AvgGowers_fromWL2.png", width = 24, height = 8, units = "in")

Subtraction Distances

wl2_garden_clim <- all_clim %>% 
  filter(parent.pop=="WL2_Garden", year>=2023) %>%
  rename_with(.cols = c(8:22), ~paste0(.x, '_WL2')) %>% 
  select(Season, WL2_year=year, cwd_WL2:ppt_coldest_month_WL2)

sub_dist <- all_clim %>% 
  filter(parent.pop!="WL2_Garden") %>% 
  group_by(parent.pop, elevation.group, elev_m, timeframe, Season) %>% 
  summarise_at(c("cwd", "pck", "ann_tmean", "ann_ppt"),
               c(mean), na.rm = TRUE) %>% #get 30 year averages
  left_join(wl2_garden_clim) %>% 
  mutate(cwd_dist=cwd - cwd_WL2,
         pck_dist=pck - pck_WL2,
         ann_tmean_dist=ann_tmean - ann_tmean_WL2,
         ann_ppt_dist=ann_ppt - ann_ppt_WL2) 
## Joining with `by = join_by(Season)`
## Warning in left_join(., wl2_garden_clim): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 3 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
sub_dist_2023 <- sub_dist %>% filter(WL2_year==2023) %>% select(parent.pop:Season, WL2_year, cwd_dist:ann_ppt_dist)
sub_dist_2024 <- sub_dist %>% filter(WL2_year==2024) %>% select(parent.pop:Season, WL2_year, cwd_dist:ann_ppt_dist)

Plot the results

#2023
sub_dist_2023 %>% 
  ggplot(aes(x=fct_reorder(parent.pop, cwd_dist), y=cwd_dist, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="CWD Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe+Season)

sub_dist_2023 %>% 
  ggplot(aes(x=fct_reorder(parent.pop, pck_dist), y=pck_dist, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="PCK Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe+Season)

sub_dist_2023 %>% 
  ggplot(aes(x=fct_reorder(parent.pop, ann_tmean_dist), y=ann_tmean_dist, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="Temp Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe+Season)

sub_dist_2023 %>% 
  ggplot(aes(x=fct_reorder(parent.pop, ann_ppt_dist), y=ann_ppt_dist, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="PPT Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe+Season)

#2024
sub_dist_2024 %>% 
  ggplot(aes(x=fct_reorder(parent.pop, cwd_dist), y=cwd_dist, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="CWD Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe+Season)

sub_dist_2024 %>% 
  ggplot(aes(x=fct_reorder(parent.pop, pck_dist), y=pck_dist, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="PCK Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe+Season)

sub_dist_2024 %>% 
  ggplot(aes(x=fct_reorder(parent.pop, ann_tmean_dist), y=ann_tmean_dist, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="Temp Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe+Season)

sub_dist_2024 %>% 
  ggplot(aes(x=fct_reorder(parent.pop, ann_ppt_dist), y=ann_ppt_dist, group=parent.pop, fill=elev_m)) +
  geom_col(width = 0.7,position = position_dodge(0.75)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_gradient(low = "#F5A540", high = "#0043F0") +
  labs(y="PPT Distance \n from WL2", fill="Elevation (m)", x="Population") +
  theme_classic() +
  theme(text=element_text(size=25), axis.text.x = element_text(angle = 45,  hjust = 1)) +
  facet_wrap(~timeframe+Season)

Merge all distances

#PREP
all_gowers_prep <- all_gowers %>% 
  select(parent.pop:Gowers_Dist) %>% 
  mutate(Season=if_else(Season=="Water Year", "Wtr_Year", "GrwSsn")) %>% 
  pivot_wider(names_from = timeframe:Year, values_from = Gowers_Dist, names_prefix = "GD_") 

wl2_gowers_avg_prep <- wl2_gowers_avg %>% 
  pivot_wider(names_from = timeframe:Season, values_from = Avg_Gowers, names_prefix = "AvgGD_") 

temp_dist_prep <- sub_dist %>% 
  select(parent.pop:Season, WL2_year, ann_tmean_dist) %>% 
  mutate(timeframe=if_else(timeframe=="historic", "Historic", "Recent"),
         Season=if_else(Season=="Growth Season", "GrwSsn", "Wtr_Year")) %>% 
  pivot_wider(names_from = timeframe:WL2_year, values_from = ann_tmean_dist, names_prefix = "TempDist_")

ppt_dist_prep <- sub_dist %>% 
  select(parent.pop:Season, WL2_year, ann_ppt_dist) %>% 
  mutate(timeframe=if_else(timeframe=="historic", "Historic", "Recent"),
         Season=if_else(Season=="Growth Season", "GrwSsn", "Wtr_Year")) %>% 
  pivot_wider(names_from = timeframe:WL2_year, values_from = ann_ppt_dist, names_prefix = "PPTDist_")

#need to add average sub dist 

#Merge
clim_geo_dist <- all_gowers_prep %>% 
  left_join(wl2_gowers_avg_prep) %>% 
  left_join(temp_dist_prep) %>% 
  left_join(ppt_dist_prep) %>% 
  mutate(WL2_Lat=38.82599, WL2_Long=-120.2509, WL2_Elev=2020.1158) %>% 
  mutate(Geographic_Dist=distHaversine(cbind(WL2_Long, WL2_Lat), cbind(Long, Lat)),
         Elev_Dist=elev_m-WL2_Elev) #Add Elevation and Geographic Distance 
## Joining with `by = join_by(parent.pop, elevation.group, elev_m, Lat, Long)`
## Joining with `by = join_by(parent.pop, elevation.group, elev_m)`
## Joining with `by = join_by(parent.pop, elevation.group, elev_m)`

Compare Growth Season and Water Year